
R version 4.5.0 (2025-04-11) -- "How About a Twenty-Six"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

- Project '~/repos/wong_bjps_2025' loaded. [renv 1.1.4]
> ## This file assesses the matches with whites using only those who were shown their DA
> ## but allowing matches across different DAs within the design
> 
> options(digits = 4, width = 132)
> 
> library(here)
> library(tidyverse)
> library(RItools)
> 
> source(here::here("Design", "nonbimatchingfunctions.R"))
> 
> load(here::here("Design", "matches_DA_new.rda"), verbose = TRUE)
Loading objects:
  matches_DA_new
> load(here::here("Data", "wrkdat_DA0_new.rda"), verbose = TRUE)
Loading objects:
  wrkdat_DA0_new
> 
> wdat0 <- wrkdat_DA0_new %>%
+   dplyr::select(-geometry) %>%
+   filter(!is.na(vm_change) & !is.na(social.capital01) & !is.na(community.resp01) & !is.na(vm)) %>%
+   droplevels()
> dim(wdat0)
[1] 650  41
> 
> wdat0 <- left_join(wdat0, matches_DA_new)
> table(is.na(wdat0$pair))

FALSE  TRUE 
  448   202 
> 
> ## Easier to think about "higher perceiver" versus "lower perceiver" when it comes to asking whether the matching worked.
> ## Further simplify the working file
> wdat1 <- wdat0 %>%
+   filter(!is.na(pair)) %>%
+   group_by(pair) %>%
+   mutate(perc_more = rank(vm) - 1) %>%
+   ungroup() %>%
+   mutate(pairF = factor(pair)) %>%
+   select(
+     perc_more, da_prop_vm_20pct_06, vm_change, pair, pairF, social.capital01, community.resp01,
+     vm, vcid, dauid, csd_pop_06, csd_prop_vm_20pct_06,
+     da_pop_dens_06, csd_urban_06, csd_pop_dens_06, community_area_km, da_pop_06
+   ) %>%
+   droplevels()
> 
> xb1_ranked0 <- xBalance(perc_more ~ da_prop_vm_20pct_06 + vm_change,
+   strata = list(unstr = NULL, pair = ~pair), report = "all", data = wdat1
+ )
> xb1_ranked0$overall
      chisquare df p.value
unstr 8.852e-05  2  1.0000
pair  1.355e-01  2  0.9345
> xb1_ranked0$results
, , strata = unstr

                     stat
vars                   Control Treatment   adj.diff adj.diff.null.sd   std.diff         z      p
  da_prop_vm_20pct_06  0.12193   0.12181 -1.153e-04         0.012777 -0.0008521 -0.009028 0.9928
  vm_change           -0.04453  -0.04451  1.811e-05         0.006998  0.0002442  0.002587 0.9979

, , strata = pair

                     stat
vars                   Control Treatment   adj.diff adj.diff.null.sd   std.diff        z      p
  da_prop_vm_20pct_06  0.12193   0.12181 -1.153e-04        0.0003145 -0.0008521 -0.36677 0.7138
  vm_change           -0.04453  -0.04451  1.811e-05        0.0011143  0.0002442  0.01625 0.9870

attr(,"originals")
[1] "da_prop_vm_20pct_06" "vm_change"          
> 
> xb1_ranked <- balanceTest(perc_more ~ da_prop_vm_20pct_06 + vm_change + strata(pair), data = wdat1, p.adjust = "none")
> xb1_ranked$overall
     chisquare df p.value
pair 1.355e-01  2  0.9345
--   8.852e-05  2  1.0000
> xb1_ranked$results[, c("Control", "Treatment", "adj.diff", "std.diff", "p"), ] # , "pair"]
, , strata = pair

                     stat
vars                   Control Treatment   adj.diff   std.diff      p
  da_prop_vm_20pct_06  0.12193   0.12181 -1.153e-04 -0.0008521 0.7138
  vm_change           -0.04453  -0.04451  1.811e-05  0.0002442 0.9870

, , strata = --

                     stat
vars                   Control Treatment   adj.diff   std.diff      p
  da_prop_vm_20pct_06  0.12193   0.12181 -1.153e-04 -0.0008521 0.9928
  vm_change           -0.04453  -0.04451  1.811e-05  0.0002442 0.9979

> 
> 
> ## An alternative for the overall test which should be fairly close to it
> library(formula.tools)
> library(coin)
> coin_fmla <- da_prop_vm_20pct_06 + vm_change ~  vm | pairF
> coin_test <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic")
> coin_test_perm <- independence_test(coin_fmla, data = wdat1, teststat = "quadratic", distribution = approximate(nresample = 1000))
> coin::pvalue(coin_test)
[1] 0.763
> coin::pvalue(coin_test_perm)
[1] 0.781
99 percent confidence interval:
 0.7455 0.8138 

> 
> coin_fmla_rank <- da_prop_vm_20pct_06 + vm_change ~ perc_more | pairF
> coin_test_rank <- independence_test(coin_fmla_rank, data = wdat1, teststat = "quadratic")
> coin_test_rank_perm <- independence_test(coin_fmla_rank,
+   data = wdat1,
+   teststat = "quadratic", distribution = approximate(nresample = 1000)
+ )
> coin::pvalue(coin_test_rank)
[1] 0.9345
> coin::pvalue(coin_test_rank_perm)
[1] 0.934
99 percent confidence interval:
 0.9111 0.9526 

> 
> pairdiffs <- wdat1 %>%
+   group_by(pair) %>%
+   summarize(
+     da_prop_vm_20pct_06 = abs(diff(da_prop_vm_20pct_06)),
+     vm_change = abs(diff(vm_change)),
+     vm = abs(diff(vm)),
+     da_pop_06 = abs(diff(da_pop_06)),
+     csd_pop_06 = abs(diff(csd_pop_06)),
+     csd_pop_dens_06 = abs(diff(csd_pop_dens_06)),
+     csd_prop_vm_20pct_06 = abs(diff(csd_prop_vm_20pct_06)),
+     # da_pop_dens_06 = abs(diff(da_pop_dens_06)),
+     community.area.km = abs(diff(community_area_km))
+   )
> 
> sapply(pairdiffs, summary)
$pair
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0    56.8   112.5   112.5   168.2   224.0 

$da_prop_vm_20pct_06
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000000 0.000266 0.001435 0.002992 0.004064 0.016267 

$vm_change
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00727 0.01448 0.01422 0.02114 0.02986 

$vm
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.010   0.090   0.170   0.247   0.310   4.190 

$da_pop_06
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    0.0    54.2   182.0   523.8   459.0 11429.0       2 

$csd_pop_06
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0   26837  257821  590654  784814 2499379 

$csd_pop_dens_06
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    0.0    35.5   282.2   846.2  1030.2  4537.6 

$csd_prop_vm_20pct_06
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00000 0.00879 0.04314 0.10093 0.17203 0.49929 

$community.area.km
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
1.20e-02 6.42e+00 4.61e+01 5.15e+02 2.43e+02 2.92e+04 

> 
> sapply(pairdiffs, function(x) {
+   quantile(x, sort(unique(c(seq(0, 1, .1), seq(.9, 1, .01)))), na.rm = TRUE)
+ })
      pair da_prop_vm_20pct_06 vm_change     vm da_pop_06 csd_pop_06 csd_pop_dens_06 csd_prop_vm_20pct_06 community.area.km
0%     1.0           0.000e+00  0.000000 0.0100       0.0          0            0.00             0.000000         1.243e-02
10%   23.3           0.000e+00  0.001581 0.0300      23.1          0            0.00             0.000000         1.333e+00
20%   45.6           9.896e-05  0.005265 0.0660      46.0      11500           24.10             0.007092         4.907e+00
30%   67.9           4.489e-04  0.008501 0.1000      84.9      55164           56.39             0.014891         8.514e+00
40%   90.2           9.154e-04  0.011661 0.1300     128.4     101246          166.82             0.031655         2.279e+01
50%  112.5           1.435e-03  0.014479 0.1700     182.0     257821          282.19             0.043136         4.612e+01
60%  134.8           2.435e-03  0.017898 0.2100     249.2     413837          399.00             0.072808         9.816e+01
70%  157.1           3.572e-03  0.019897 0.2710     375.5     685513          858.66             0.130674         1.783e+02
80%  179.4           5.411e-03  0.022633 0.3400     586.4    1142080         2074.33             0.191346         3.836e+02
90%  201.7           9.019e-03  0.026221 0.4800    1122.8    1925240         2904.36             0.302753         8.628e+02
91%  203.9           9.086e-03  0.026390 0.4993    1324.0    1925240         2912.35             0.306581         1.177e+03
92%  206.2           9.411e-03  0.026802 0.5100    1454.0    2016255         2931.77             0.328060         1.215e+03
93%  208.4           9.916e-03  0.027148 0.5278    1565.4    2124912         3090.47             0.331926         1.530e+03
94%  210.6           1.079e-02  0.027474 0.5824    1804.7    2235158         3118.48             0.336357         1.740e+03
95%  212.9           1.113e-02  0.028284 0.6740    2072.7    2300202         3193.78             0.369550         1.949e+03
96%  215.1           1.123e-02  0.028435 0.7632    2452.2    2390318         3303.52             0.377446         2.049e+03
97%  217.3           1.163e-02  0.029130 0.9117    2638.6    2428464         3393.22             0.399906         2.658e+03
98%  219.5           1.268e-02  0.029422 1.0362    3827.6    2444732         3464.38             0.416627         3.909e+03
99%  221.8           1.409e-02  0.029688 1.3231    6011.4    2479360         3902.14             0.433720         4.999e+03
100% 224.0           1.627e-02  0.029860 4.1900   11429.0    2499379         4537.65             0.499291         2.923e+04
> 
> ### facts for the body.tex
> 
> # number of people in final design
> nrow(wdat1)
[1] 448
> # number of DAs
> length(unique(wdat1$dauid))
[1] 443
> 
> system("touch Design/match_assess_DA_new.done")
> 
